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The deformation of an initially spherical capsule, freely suspended in simple shear flow, can 
be computed analytically in the limit of small deformations [D. Barthes-Biesel, J. M. Rallison, The 
Time-Dependent Deformation of a Capsule Freely Suspended in a Linear Shear Flow, J. Fluid Mech. 
113 (1981) 251-267]. Those analytic approximations are used to study the influence of the mesh 
tessellation method, the spatial resolution, and the discrete delta function of the immersed boundary 
method on the numerical results obtained by a coupled immersed boundary lattice Boltzmann finite 
element method. For the description of the capsule membrane, a finite element method and the 
Skalak constitutive model [R. Skalak et al., Strain Energy Function of Red Blood Cell Membranes, 
Biophys. J. 13 (1973) 245-264] have been employed. Our primary goal is the investigation of the 
presented model for small resolutions to provide a sound basis for efficient but accurate simulations 
of multiple deformable particles immersed in a fluid. We come to the conclusion that details of the 
membrane mesh, as tessellation method and resolution, play only a minor role. The hydrodynamic 
resolution, i.e., the width of the discrete delta function, can significantly influence the accuracy 
of the simulations. The discretization of the delta function introduces an artificial length scale, 
which effectively changes the radius and the deformability of the capsule. We discuss possibilities 
of reducing the computing time of simulations of deformable objects immersed in a fluid while 
maintaining high accuracy. 

PACS numbers: 47.ll.-j, 47.57.-s, 47.63.-b 



I. INTRODUCTION 

Understanding the hydrodynamics of blood is certainly 
one of the major motivations for the simulation of de- 
formable particles immersed in a fluid. Aside from the 
desire to investigate this aspect of fundamental research 
in more detail, there are many relevant applications in 
biology and medical sciences. The ultimate goal is the 
simulation of the human microcirculation up to the cen- 
timeter scale including the full dynamics of the cells, their 
interactions with each other and the blood vessel walls, 
and the impact of their microscopic properties on the 
macroscopic behavior of blood over the entire shear rate 
range. The complexity and scale-bridging of the coupled 
system of hydrodynamics and cell membrane dynamics 
requires numerical approaches to obtain meaningful re- 
sults. 

Capsules are elastic membranes filled with a fluid. The 
investigation of the dynamical response of capsules in flu- 
ids is not trivial since the hydrodynamic properties of 
the internal and external fluids, the constitutive model 
of the capsule membrane, and its shape affect the out- 
come. Fortunately, there are analytic solutions available 
for small deformations of single capsules in shear flow 
P, Q • Experimental results for the behavior of artificial 
capsules have been obtained by Chang and Olbricht Q 
and Walter et al. [||. Pozrikidis Q numerically stud- 
ied those problems using the boundary element method 
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(BEM). This method, in different formulations, has also 
been employed by, e.g., Kraus et al. @, Ramanujan and 
Pozrikidis Q, Diaz et al. Q, Barthes-Biesel et al. 
and Lac et al. [Io| • The BEM is only valid in Stokes flow 
and cannot be applied to flow situations with inertia. 
Fluid-capsule interactions using the immersed boun dary 
method have been studied by Egglcton and Popel [ill ]. 
Sui et al. [l2j have implemented a combined immersed 
boundary lattice Boltzmann method with grid refinement 
to study the transient deformation of capsules in simple 
shear flow at high resolutions. Pozrikidis [l3| consid- 
ered bending resistance in capsule simulations. He states 
the importance of the presence of bending stiffness in bi- 
ological cells. Simulations of multiple blood cells have 
been performed by, e.g., Dupin et al. (l4| and Doddi and 
Bagchi pi. 

Although the processing power of present-day comput- 
ers is large compared to that one or two decades ago, 
computing resources arc still limited. In simulations of 
multiple deformable particles in flow, the spatial resolu- 
tion must be kept sufficiently small in order to limit the 
computing time to a reasonable period. This calls for effi- 
cient numerical methods which are capable of capturing 
the dominant physical properties of the problem, even 
at small resolutions. Therefore, we aim at a better un- 
derstanding of the behavior of the simulation technique 
employed in this paper, especially at smaller resolutions. 

The lattice Boltzmann method (LBM) is a compara- 
bly new method to solve the full Navier-Stokes equations 
pL6l - l20j . The starting point is the lattice Boltzmann equa- 
tion which is an approximate and discrctized form of the 
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Boltzmann equation. Virtual particles, also called popu- 
lations, arc moving on a lattice and collide at the fixed 
lattice nodes of the regular grid. In the macroscopic 
limit, the Navier- Stokes equations can be recovered from 
the well-defined collision rules of those particles. This 
method is particularly straightforward to implement, and 
it has proven to be accurate and applicable to many hy- 
drodynamic problems, e.g., f2ll - |24j ] . 

Peskin developed the immersed boundary method 
(IBM) [H [IH to model blood flow in the heart. The 
strength of this method is that the immersed material is 
intrinsically deformablc and that the Navier-Stokes equa- 
tions do not need to be modified in the presence of the 
material, except for the inclusion of a body force. The 
IBM has especially attracted the attention of the scien- 
tists due to its capability to model thin and elastic mem- 
branes and therefore red blood cells and capsules in ar- 
bitrary external flow fields. In an attempt to benchmark 
and test the accuracy and convergence behavior of IBM, 
the IBM has been applied to stiff objects [27l - l30| . How- 
ever, there arc numerical limits to the applicability of 
IBM to stiff materials, and it requires some effort to do 
so. On the contrary, the lattice Boltzmann bounce-back 
(BB) scheme unrolls its benefits in the case of stiff obsta- 
cles in flow [3 I30l - l34j , and it is much more demanding 
to capture deformable objects with BB. 

In order to efficiently simulate deformable membranes 
immersed in a fluid, we use an approach combining the 
LBM as fluid solver, the IBM for the coupling of the 
fluid and the membranes, and a finite element method 
(FEM) for the computation of the membrane response to 
deformations (IBLBFEM). This approach has been em- 

ted successfully by other scientists, e.g., Zhang et al. 
or Sui et al. [12|. A similar method, with another 
Navier-Stokes solver, has also been applied by Eggleton 
and Popel [TTj] and Doddi and Bagchi The advan- 
tage of this combined IBLBFEM is that the computa- 
tions of the fluid and membranes are decoupled and that 
the meshes of the fluid and the membranes do not have 
to match. No remeshing is required for the membranes. 
The implementation is straightforward, and the method 
is powerful to simulate 0(100) deformable particles in 
flow solving the full Navier-Stokes equations and obtain- 
ing velocity, pressure, and shear stress information locally 
and at finite Reynolds numbers. 

We use an explicit IBM coupling scheme which is ef- 
ficient in terms of computing time. However, the no- 
slip condition at the membrane surface is not perfectly 
obe yed , and a drift of the capsules' volume can occur 
[36l |37(. For the simulation of stiff objects in flow (which 
is not the case here) , modified IBM schemes exactly obey- 
ing the no-slip condition are known [H, |39[ . It is also of- 
ten argued that the IBM is unstable in the limit of high 
stiffness [4(| ■ There are methods to increase stability us- 
ing either implicit or semi-implicit methods [37l l4lj . Due 
to the softness of the capsules and the relatively short 
simulations in the present paper, instabilities do not oc- 
cur, and the volume drift is negligible. We do not in- 



tent to comment on IBM-related stability issues or an 
improved implementation of the no-slip condition in this 
work. 

When it comes to the discretization of the capsule, 
the question arises whether a structured or unstructured 
mesh should be used and how this mesh should be cre- 
ated. Structured meshes usually contain coordinate sin- 
gularities whereas on unstructured meshes gradients have 
to be approximated. Diaz et al. @ and Lac et al. 1(| have 
used structured meshes whereas Kraus et al. [6|, Navot 
[43 | . and Ramanujan and Pozrikidis Q have employed 
unstructured grids for the capsule. We use an unstruc- 
tured mesh tessellation in this paper. 

This paper is not targeted at finding new physics of 
capsules in shear flow. We rather wish to better under- 
stand the behavior of the combined IBLBFEM and the 
membrane tessellation on the accuracy and the numeri- 
cal efficiency of the simulation of deformable capsules at 
small resolutions. This is the main difference between 
our recent effort and the work by Sui et al. [l2| who have 
not discussed the behavior of their numerical method at 
small resolutions. Of special interest are the impact of the 
IBM interpolation stencil, the mesh tessellation, and the 
ratio between the average mesh node distance and the 
lattice constant of the LBM grid. We will also shortly 
discuss the importance of a well-chosen BGK LBM re- 
laxation parameter. Those considerations are important 
for a correct setup of efficient simulations containing a 
large number of deformable objects with relatively coarse 
meshes. Therefore, our investigations are hoped to be 
useful for future simulations of multiple deformable par- 
ticles in flow employing the IBM. We simulate the time 
evolution of the deformation of an initially spherical cap- 
sule freely suspended in an unbound simple shear flow. 
The interior and exterior fluids have the same properties, 
and they are Newtonian. The capsules have no bending 
resistance. For comparison, the approximated analytic 
steady-state solutions by Barthes-Biesel [l[ and Barthes- 
Biesel and Rallison 0] are used. 

The LBM will be shortly presented in Sec. Ill At fol- 
lowed by an overview of the membrane model and the 
used F EM in Sec. iHBl The IBM is briefly presented 
in Sec. Ill C( and the mesh influence is discussed in Sec. 
Ill Dl In Scc lIII) the theory of small capsule deformations 
is shortly outlined. The simulations, results, and discus- 
sions can be found in Sec. lIVI followed by the conclusions 
in Sec. EJ 



II. NUMERICAL METHODS 

The simulation algorithm consists of three major com- 
ponents: the fluid solver, the membrane model, and the 
coupling of fluid and membrane. For the fluid solver, 
we have used the D3Q19 Bhatnagar-Gross-Krook (BGK) 
lattice Boltzmann method (LBM). The membrane dy- 
namics is derived from a constitutive model, and the 
strains in the capsule are evaluated using a finite ele- 
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mcnt method (FEM). The interaction of the fluid and 
the membrane is captured by the immersed boundary 
method (IBM). 

The overview of LBM is presented in Sec. Ill Al the 
membrane model is outlined in Sec. Ill Bt and IBM is 
shortly covered in Sec. Ill CI Since also the discrete par- 
ticle mesh plays a role in the simulations, its properties 
are briefly discussed in Sec. Ill Dl 

A. Lattice Boltzmann method 

In the last two decades, the LBM has become a com- 
petitive Navier-Stokes solver with increasing prominence 
among scientists in the field of computational fluid dy- 
namics fl6l - [l8l . HH, HiJ . Its strength is based on its sim- 



ple coding and, since LBM is an automaton, its locality, 
making it intrinsically parallelizable. 

The lattice Boltzmann equation (LBE), Eq. ([1]), can be 
regarded as a discretized form of the Boltzmann equation. 
On the other hanchit is an extension of the lattice gas cel- 
lular automaton [45| . In contrast to conventional Navier- 
Stokes solvers, the LBE is not the discretized form of the 
Navier-Stokes equations. While conventional methods di- 
rectly solve the Navier-Stokes equations in terms of the 
pressure p and the velocity u, the LBM introduces a num- 
ber of q populations /, (i = 0, . . . , q — 1) streaming along 
a regular lattice (lattice constant Ax) in discrete time 
steps. Those populations can be regarded as mesoscopic 
particle packets propagating and colliding. 

The evolution of the populations fa is given by the 
LBE, which takes the form 



fi{x + CiAt, t + At) - fi(x, t) = — (fa{x, t) - fP(x, t)) + FiAt 

T 



(1) 



in the BGK approximation. The dimensionless relax- 
ation parameter t of the fluid is connected to the speed 
of sound c s and the kinematic viscosity v by v — c 2 s (r — 
1/2) At, where c s = Ax/ At holds. At each time 

step t, the populations propagate along the q discretized 
velocity vectors Cj to the next neighbors. At those points, 
they collide according to the right-hand side of Eq. (J} . 
The significance of Fi is explained below. The equilib- 
rium populations are given by 

/f 5 = Wi p ^1 + 3cj -u + ^(ci- u) 2 - |it • v^j . (2) 

This is closely related to the truncated form of the 
Maxwell distribution which is a very good approximation 
for small Mach numbers. The q factors to, are the lattice 
weights, depending on the underlying lattice structure. 
Their choice ensures the isotropy of the fluid, a neces- 
sity to solve the Navier-Stokes equations asymptotically. 
In the present paper, we use a 3D model with 19 ve- 
locities, designated D3Q19. The lattice structure, the 
corresponding velocities Ci, and the lattice weights Wi 
are introduced in [l||. A sketch of the D3Q19 lattice is 
shown in Fig. Q] 

A body force density / can be incorporated via Fi in 

Eq. © EMEU, 




This force density is particularly important for the cou- 
pling of the fluid and the immersed membranes, but it is 
also commonly used to include gravity. More details are 
given in Sec. Ill CI 

Finally, the macroscopic properties of the fluid have to 
be extracted from the populations fi. The density and 



the velocity can directly be recovered by computing the 
zcroth and first moments: 

= (4) 

i 

x , At 
pu = 2_^Cifi + —f, (5) 

i 

at each fluid lattice node. The deviatoric shear stress ten- 
sor <j can also be computed from the populations locally 

The simple shear flow required for the present simula- 
tions can be realized by using the bounce-back method 
for moving walls [l8| . The fluid is fully periodic along the 
x- and y-axes (velocity and vorticity directions, respec- 
tively), but it is bound by two plane walls at z = ±H/2, 
where H is the distance of the walls. Moving the walls 
in the x-direction with velocity ±m„ (u w > for z > 0), 
the shear rate of the fluid is 7 = 2u w /H. 

In order to obtain accurate predictions, the limits of 
the validity of the LBM must be acknowledged. The 
slip velocities u w have to be chosen sufficiently small so 
that the LBM operates in the small Mach number limit. 
An even more stringent restriction of the wall speed is 
given by the Reynolds number. The theory of small de- 
formations of a capsule in simple shear flow is only valid 
in Stokes flow, Re < 1. The correct choice of the re- 
laxation parameter r also influences the accuracy of the 
simulations 0, 5H , especially in combination with the 
IBM ■ The LBE intrinsically contains the partial time 
derivative du/dt of the Navier-Stokes equations. Hence, 
Stokes flow can only be reached asymptotically. 

A more detailed presentation of the LBM can be found 
in the literature, e. g., i n the monographs by Succi [l!| or 
Sukop and Thorne [20j]. 
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Figure 1. Sketch of the D3Q19 lattice structure. All velocity vectors are located in at least one of the three coordinate planes 
(light gray). The velocity vectors either point to the next neighbors along the coordinate axes (ci_6, black arrows) or to the 
next but one neighbors (c7_is, dark-gray arrows). The zero velocity co is not shown. 



B. Membrane model and force computation 

The IBM algorithm requires knowledge of the forces 
at the nodes of the tessellated membrane (cf. Sec. Ill CI) . 
For a hyperelastic material, i.e., negligible viscous and 
plastic forces, the shear forces can be computed from a 
constitutive model for the areal strain energy density w . 
The contributions to the total energy W can generally 
be written in the form W = W s + W B + W A + W v , 
where the superscripts denote strain, bending, surface, 
and volume contributions, respectively. 

The areal strain energy density w obeying W s — 
J dAw s (dA is the surface element) can only depend 
on the invariants I\ = Af + \\ — 2 and I2 = A^A| — 1 for 
a thin membrane with isotropic and homogeneous clastic 
properties, i.e., per definition w s is invariant under rota- 
tions and translations. Ai and A2 are the local principal 
in-plane stretch ratios. Deformations of biological cells 
can be large, and thus the linear strain-stress approxi- 
mation is not justified in general. Skalak et al. [48| have 
suggested an energy model which is able to reproduce 
experimental data of red blood cells at both small and 
large strains, 
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The surface elastic shear modulus k s and area dilation 
modulus k a control the strength of the membrane re- 
sponse to deformation (shear and dilation). A commonly 
used model is the neo-Hookean law which is equivalent 
to the zero-thickness shell membrane proposed by Ra- 
manujan and Pozrikidis for small deformations. An- 
other constitutive model has been proposed by Navot 
1421 . More information about those constitutive laws can 



be found in the literature, e.g., HdH and will not be dis- 
cussed here. For all our simulations, we have employed 
the Skalak membrane model, Eq. (|6|). 

The bending energy W B can have local and non-local 
contributions. We have not considered any bending en- 
ergy in the present simulations, i.e., W B = since it is 
not considered in the analytic investigation by Barthes- 
Biesel and Rallison 0. However, it has been thor- 
oughly discussed in the literature that a bending resis- 
tance has to be included whenever strong local curva- 
tures appear. Elsewise, the membranes can buckle or 
collapse @, El, ED, This is especially the case for 

more complex geometries and strong deformations as in 
the case of red blood cells. More details about the form of 
the bending energy are provided in Canham , Helfrich 
[5lT ] , Svetina and Zeks [52j , and Gompper and Schick [HI] . 

The total volume and surface of the membrane may 
be restricted. This can be formulated by defining global 
volume and surface energies, W v and W A . Those ener- 
gies are minimum if the volume and surface equal their 
corresponding equilibrium values [53j . In the present sim- 
ulations, we have neither employed a volume nor a sur- 
face energy, i.e., W v = W A = 0. Although the IBM is 
not perfectly volume-conserving, cf. Sec. Ill Cl the volume 
drift is almost negligible in most of the present simula- 
tions of a single capsule. The reason is the relatively short 
duration of the simulations. We stress that volume and 
surface energies may have to be considered to improve 
numerical stability in longer simulations with complex 
flow fields, large local shear rates, and coarse mesh reso- 
lutions. 

The capsule membrane is numerically described by a 
number Nf of flat triangular face elements, which remain 
flat even at large deformations. While the deformation 
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Figure 2. Illustration of (a) the equilibrium face (denned by lo, 1' , and tpo), (b) its deformed shape (accordingly defined by I, 
I' , and ip), and (c) both transformed to the same xy- plane. The displacement vector v\ is identically zero, and the other two 
are shown in subfigure (c). The deformation state (Ai, A2) of the face is then uniquely defined. 



state is a property of the faces, the membrane forces have 
to be known at the corners of the faces (nodes) . The first 
step in the computation of the strains A1.2 of a given 
face element is the identification of the displacements Vi 
(i = 1, 2, 3) of the nodes. The deformed and undeformed 
elements are transformed to a common plane (here: xy- 
plane) in such a way that the edges 1' and I' are aligned, 
cf. Fig. [5J There is no restriction since translations and 
rotations do not change the energy of an element. The 
basic assumption is that the displacement gradient ten- 
sor {D a p) w (5 a j3 + dpv a ) is spatially constant over the 
entire face element. This can be realized by introduc- 
ing a linear shape function Ni{x,y) = a^x + biy + c.j 
(i = 1,2,3) for each node. The coefficients are found 
by letting Ni(xj 7 yj) — Sij = 1,2,3), i.e., each shape 
function Ni is unity at the location of the correspond- 
ing node i, but zero at the two nodes other than i. The 
linear displacement field of the face element can then be 
written as 



from the equations 



v(x, y) = N1V1 + N 2 v 2 + N 3 v 3 , 



(7) 



and the displacement gradient tensor D can be com- 
puted. Its components do not depend on x or y, but on 
the shape function coefficients a; and bi which are fixed 
by the shape of the undeformed element only, i.e., a, and 
bi are constant in time for each node in the face. It is 
straightforward to show that the displaceme nt g radient 
tensor then has the form {D a p) = (g h c ) with |53j 



A X A 2 
\\ + \ 2 2 = a 2 



2 2 
a c , 



(11) 
(12) 



since \\ + X 2 . = tr D T D and XfX^ = det D T D. Note that 
the product D T D is rotationally invariant. For more 
details, we refer to Charrier et al. 1541. Shrivastava and 
Tang [55| . and Gompper and Schick [53| . 

The total energy in the present simulations has shear 
contributions only, W = W s . The strain energy W s is 
computed from the areal energy density w s , Eq. ([6]), and 
the local reference area A of the membrane face elements 

MM, 



faces 



(13) 



Once the energy of the membrane is known, the forces 
acting on the fluid exerted by node i at position Xi can 
be computed from the principle of virtual work, 



dW{xj) 
dxi 



(14) 



This procedure is equivalent to the approach explained 
in details in [541151. 



/ 

a= r , 
to 
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I' sin tfi 
1' sin ip ' 



COS IfiQ 



(8) 
(9) 
(10) 



Here, I and I' are the lengths of two arbitrary edges of the 
face, and ip is the angle between those edges, cf. Fig. [5] 
The zero index (Iq. 1' , and ipo) denotes the undeformed 
values. The current deformation of a face is evaluated 



C. Immersed boundary method 

The IBM was originally proposed by Peskin [2^, [2(| . 
The basic idea is to couple the Eulcrian coordinate sys- 
tem of the fluid lattice and the arbitrary Lagrangian co- 
ordinate system of a surface which is not conform to the 
regular lattice. The IBM is a front-tracking coupling 
method. 

Caused by its deformations, the membrane exerts a 
force Fi(t) on the fluid at time step t. The fluid lattice 
nodes have fixed positions X and the membrane nodes 
i are located at Xi(t). In the discretized description, the 
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Figure 3. Illustration of the immersed boundary method in 
two dimensions. The membrane mesh (light gray) moves on 
top of the fixed fluid lattice (black). If the 2-point interpo- 
lation stencil is used, only the four lattice nodes enclosed by 
the dashed square with side length 2Ax are required for the 
spreading and interpolation steps of membrane node i located 
at Xi(t) (dark-gray). 



Eulcrian body force density f(X,t) is computed from 
the Lagrangian force Fi (t) by the spreading operation 



f{x,t)=Y,m)s{x-x i {t))- 



(15) 



The lattice force density f(X,t) is then used in the lat- 
tice Boltzmann equation, Eq. ([1]), via the force coupling, 
Eq. ([3]). The kernel S(X — Xi{t)) is a discretized Dirac 
delta function with a finite support. Peskin (2(| has 
shown that this function has to obey some basic prop- 
erties to maintain momentum and angular momentum 
conservation. Still, the width of the function is not re- 
stricted a priori, and it can be considered as a free pa- 
rameter of the IBM. We use the common decomposi- 
tion S(r) = 4>{x)4>(y)(f)(z). Among others [26], (5(|, the 
most popular interpolation functions <p n with a support 
of n = 2, 3, 4 lattice nodes along each coordinate axis are 



h(r) = 
h(r) = 



Mr) 




(l + VI - 3r 2 ) 

h - 3|r| - y/-2 + 6\r\ - 3r 2 



3 - 2\r\ + v/l + 4|r|-4r 2 



5 - 2\r\ - ^-7+12^1 -4r 2 



o<H<| 

k<\r\<h 
1<H 

< \r\ < 1 

1 < \r\ < 2 

2 < Irl 



(16) 
(17) 

(18) 



Completing the coupling of the fluid and the membrane, 
the new velocities u^t+l) of the membrane nodes i are 
computed in the interpolation step, using the new lattice 
velocities but the old node positions, 

Ui(t + At)=Y^ u(X, t + At) S(X - Xi(t)). (19) 
x 

Here, the no-slip condition is assumed to be valid at the 
location of the membrane, i.e., the membrane moves with 
the ambient fluid velocity. The interpolation functions in 
Eqs. (fT5|) and (|T9|) are the same. The principle of the IBM 
is illustrated in Fig. [31 Finally, the membrane nodes i are 
advected explicitly by the Euler rule 

Xi(t + At) = Xi(t) + m(t + At) At (20) 

We have found that the Adams-Bashforth scheme 

Xl (t + At) = Xi {t) + (^Ui{t + At) - ^Ui(t)j At (21) 

which has been used by Doddi and Bagchi [H?} does not 
change the results for simulations of short duration. Yet, 



it provides additional accuracy for long-time simulations 
since it is a second-order scheme. 

Although the no-slip condition at the membrane sur- 
face can be exactly fulfilled in the continuous limit, this 
is not the case in the discretized, explicit version of IBM. 
Even for an incompressible velocity field, the standard 
interpolation algorithm shown in this section does not 
assure that the volume of a closed membrane remains ex- 
actly constant in time. This problem has been recognized 
early, and improved immersed boundary approaches have 
been proposed for example by Peskin and Printz [36j and 
Wu et al. Due to the comparably short simulation 
times in the present paper, there is no need to counteract 
the volume drift. 

In conclusion, each time step of the combined 
IBLBFEM scheme consists of the following sub-steps (we 
set At = 1). 

1. At the beginning of time step t, the membrane node 
positions Xi(t) and the entire fluid state u(X,t), 
p(X,t) are known. From the displacements of the 
membrane nodes, the forces -Fi(i) are computed us- 
ing the FEM (SecHH]). 



7 




n n 



Figure 4. Illustration of the mesh subdivision. The face de- 
fined by the nodes m, Tij, and nt (dark-gray) is subdivided. 
First, the edges are halved and new nodes ni, n m , n n are cre- 
ated. The six nodes are connected in such a way that four 
faces of equal area are produced (dashed lines). Finally, the 
nodes n;, n m , and n n are radially shifted, until they are lo- 
cated on the sphere enclosing the body. The four final faces 
are shown in light gray. 



2. The membrane forces Fi(t) are spread to the Eu- 
lerian grid via IBM, Eq. (fT5|) . and the body force 
density f(X,t) is obtained. 

3. f(X,t) is used in the LBM to compute the new 
state of the fluid, u(X,t + 1), p{X,t + 1) (Sec. 

mm- 

4. The new velocities Ui(t + 1) of the membrane nodes 
are computed in the framework of IBM, Eq. (fi"9l) . 

5. The new positions of the membrane nodes Xi(t + 1) 
are found by evaluating Eq. (|20p . 

6. Information on the membrane and fluid state after 
time step t may be written to the disk, using Xi {t + 
1), u(X,t + 1), etc. Get back to the first sub-step 
and recompute for time step t + 1 . 



D. Membrane tessellation 

The question arises whether the detailed properties of 
the mesh tessellation have a significant impact on the 
quality of the simulation results or not. There are dif- 
ferent approaches for the generation of a spherical mesh. 
We will focus on three of them: 

1. tessellation of an implicit surface using the CGAL 
libraries [59j . 

2. finite clement mesh generation using Gmsh [o0| . 

3. successive subdivision, starting from a coarse mesh 
of high symmetry. 



The CGAL libraries [59| allow the user to tessellate the 
surface defined by the zero level set of any implicit func- 
tion F(x, y, z). For a sphere with radius r, this function 
reads F(x,y,z) = x 2 + y 2 + z 2 — r 2 . One has control over 
the number of faces and the distance between neighboring 
nodes, but the mesh suffers from a reduced homogeneity 
and isotropy. The advantage of this method is that any 
implicit surface can be tessellated without much effort. 

Gmsh [6(| is not able to tessellate implicit surfaces, 
but the user can construct geometric objects like spheres. 
The surface of those shapes can then be tessellated. Also 
here, the high isotropy of the initial sphere is not com- 
pletely captured by the mesh. 

A method to produce a spherical mesh of high homo- 
geneity and isotropy is subdivision of a highly symmetric 
mesh of low resolution. We have used a regular icosa- 
hedron, one of the five Platonic solids. It has 20 equi- 
lateral triangles as faces, 12 nodes and 30 edges of equal 
length. The numbers of nodes N n and faces Nf of any 
closed surface consisting only of triangles are related by 
2N n = Nf + 4. Ramanujan and Pozrikidis Q and Sui 
et al. [l2j have used a similar approach, starting from 
a regular octahedron. The subdivision scheme starts at 
creating a new node at the middle of each edge. Those 
initially 30 new nodes (in case of an icosahedron) are 
then radially shifted until they are located on the cir- 
cumsphere of the body and connected to form additional 
faces. This procedure can be repeated numerous times. 
It is illustrated in Fig. 2] It has to be noted that each 
subdivision step increases the number of faces according 
to NJ 1 = N® ■ 4 m , where m is the number of subdivisions 
and N® = 20 for an icosahedron and 8 for an octahedron. 
Although the resulting mesh has surpassing properties in 
terms of edge length, face area, and angle distributions, 
one cannot create a mesh with an arbitrary number of 
faces. This restriction can somewhat be relaxed by start- 
ing from another body of lower symmetry and another 
number of faces. 

In Tab. U the properties of the meshes created by the 
presented methods are listed and compared for the case 
of a sphere with Nf 1280 faces. The meshes are il- 
lustrated in Fig. [5] Obviously, the subdivided mesh has 
the smallest scatter in face area, edge length, normal-to- 
normal angle, and edge-to-edge angle distributions. By 
default, we use the mesh obtained from the icosahedron 
if not elsewise stated. 



III. THEORY 

We assume that the ambient fluid and the fluid inside 
the capsule are Newtonian and have the same properties, 
especially the same density (p = p m = p ou t) and viscosity 
(A = 77i„/?7out = 1 and y = n„ = Vout)- This commonly 
adopted assumption [y, LL2, |42j significantly simplifies the 
computations without losing too much generality. In the 
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Table I. Properties of spherical meshes, created with CGAL, Gmsh, and a successive subdivision of a regular icosahedron. The 
spheres have Nf ~ 1280 nodes each. N„ n and N„ n are the smallest and largest number of neighboring nodes to any node found 
for the given mesh. A is the face area, I is the edge length, ip n is the angle between neighboring face normals, and ip e is the 
angle between neighboring edges (both being members of a common face). The bar indicates the average of a quantity taken 
over the entire mesh, a denotes the standard deviation. Clearly, the mesh based on the icosahedron has superior quality in 
terms of isotropy and homogeneity. 



quantity CGAL Gmsh subdivision 

W f 1278 1284 1280" 

N < N> 4 10 4 8 5 6 

o A /A 25.8% 21.2% 8.6% 

a t /I 19.0% 14.3% 6.5% 

a v Jq> n 42.5% 26.2% 15.9% 

er v J(pe 25.1% 17.2% 9.3% 



limit of Stokes flow, the capsule Reynolds number 



Re 



(22) 



is small, and inertia effects can be neglected. The only 
physical parameter left is the dimcnsionless shear rate 



G 



ka 



(23) 



where 7 is the shear rate of the unperturbed ambient 
fluid, and r and k s are the radius and the surface elastic 
shear modulus of the initially spherical capsule. 

Due to the presence of the external shear flow, the 
membrane deforms. The generated membrane tensions 
oppose the shear forces exerted by the fluid. For G <C 
1, a small membrane deformation suffices to compen- 
sate the shear forces, and the capsule shape is only 
slightly perturbed. In a simple shear flow, after an ini- 
tial transient, steady tank-treading motion of an initially 
spherical capsule, deformed to an ellipsoid, is observed 
[E HI, & S IHl ■ Tank-treading has been described first 
by Schmid-Schonbein and Wells [Hj]. The shape and dy- 
namics of the capsule can then be defined by three con- 
stant parameters: the Taylor deformation parameter D, 



the inclination angle 8, and an angular velocity uj. The 
stationary geometry is shown in Fig. [6] 

The inclination angle 6 of the membrane is taken be- 
tween its largest semiaxis r a and the x-axis, the direction 
of the fluid velocity. Dupin ct al. |l4| have extracted 9 



from an ellipsoid fit of the membrane cross-section on 
the xz-plane. A common alternative is the comparison 
of the capsule shape with an ellipsoid with the same in- 
ertia tensor I. Following Ramanujan and Pozrikidis 0, 
its components are given by 



In 



^ faces 



(24) 



where is the centroid of face i and rii its normal. The 
diagonalizcd inertia tensor of an ellipsoid with unit mass 
and constant density is / = diag(r^+r^, r^+rf , r^+r^)/5. 
r a and r c are the largest and smallest semiaxis of the el- 
lipsoid, and rb is the intermediate semiaxis. Assuming 
that the symmetry axes of the deformed capsule coincide 
with the principal axes of the inertia tensor, the ellip- 
soid's principal semiaxes and the inclination angle can 
be computed. For small deformations, this is an excel- 
lent approximation. 
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Figure 6. Sketch of the tank-treading geometry. The capsule 
cross-section is shown in the xy-plane. It is deformed with 
major and minor semiaxes r a and r c . The inclination angle 8 
is taken between the major semiaxis and the x-axis (velocity 
direction of the external flow) . The membrane rotates about 
its spatially fixed shape with angular velocity tj. The flow 
direction of the unperturbed ambient fluid is also shown (dark 
gray arrows). 



The Taylor deformation index is defined as 



D 



> 0. 



(25) 



It is a measure for the deformation of the capsule, and 
D = holds for a sphere. In the case of small deforma- 
tions in Stokes flow (G, Re <C 1), the shape of the capsule 
can be described analytically |2| . In simple shear flow in 
the steady state, the relation between the dimensionless 
shear rate G and the deformation parameter D is 



D = - 



3«2 + 4«3 



-G + 0{G 2 ). (26) 



4 a\(Za 2 + 5a 3 ) + 2a 3 (a 2 + a 3 

The coefficients a±, a 2 , and a 3 can be extracted from the 
constitutive model, Eq. (|6|), by expansion, 



w s /k s = w + «i A a + -(a 2 + ai)A^ + a 3 (A b - A«) + 0(A 3 a ,A a A b , A 2 b ), (27) 



where A a = ln(AiA 2 ) and 2A b = A? + A§ - 1. This leads 
to a.\ — 0, a 2 = 2/3, and a 3 = 1/3 for the neo-Hookcan 
law and for Skalak's law, Eq. ©, with k s = k a . Barthes- 
Biesel [l[ has also found an analytic expression for the 
deviation of the inclination angle at small G, 

6 l = ^l= l ±G (28) 

7T 7T 8 

for a neo-Hookean membrane and Skalak's law with k s = 
k a . In the limit G —> 0, the angle approaches 8 = 7r/4. 
The opposite angle 9° is the relevant angle in this problem 
since it is proportional to G. It is the deviation from the 
inclination angle of a stiff sphere, where 9 = it/ A. 

In the steady state, the membrane nodes rotate about 
the fixed shape of the capsule. This tank-treading be- 
havior can be quantified by the rotation period T or the 
angular velocity U) in the steady state. Kraus et al. @ 
have measured the time between two successive vertex 
crossings of the xy-plane. This time corresponds to half 
the rotation period. However, this value is time inte- 
grated and cannot resolve possible numerical fluctuations 
of the angular velocity. Another approach is the direct 
computation of the angular velocity of a membrane node 
by 0Ji = Atfi/At, where Acp^ is the angle swept by that 
node projected on the xz-plane during time At. For a 
stiff sphere, the angular velocity is uj = 7/2, but u> < 7/2 
holds for dcformablc capsules. 

Additionally to the approach using the inertia tensor, 
we have computed the deformed capsule shape by a linear 



I 

fit of the membrane node positions. Minimizing 

nodes 



with respect to the fit parameters Q y , £ z , and Cxz and 
diagonalizing the matrix 

/ Cx -CxA 

Q = ( y , (30) 

\-Cx, c, / 

one can directly compute the semiaxes, the deformation 
parameter, and the inclination angle about the y-axis. 
For convenience, we have not allowed rotations about 
the x- or the z-axis in Eq. (|29[) . Comparing the results 
obtained from the inertia tensor and the linear fit, we 
observe that both the inclination angle and the deforma- 
tion parameter are virtually identical. For that reason, 
we will show the data obtained from the inertia tensor 
only. However, the values of the semiaxes r a , r&, and 
r c arc slightly underestimated by the inertia tensor. The 
reason is the discretization of the mesh and the use of flat 
triangular elements. Since the semiaxes are not explic- 
itly required for the characterization of the deformation 
and the values of D and 9 are correct, we have not at- 
tempted to improve the results obtained by the inertia 
tensor method. 
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IV. SIMULATIONS AND RESULTS 

All distances are made dimensionlcss using the lattice 
constant Ax as characteristic length scale. As mentioned 
above, from the constitutive law, Eqs. (O and (|27[) . the 
expansion parameters a\ = 0, a-i = 2/3, and «3 = 1/3 
follow for k s = k a . Inserting these values in Eq. (f26|) . one 
obtains D/G = 25/4 in the linear regime. We restrict 
ourselves to the Skalak model, Eq. with k s = k a 
in all the simulations. If not otherwise stated, the lattice 
Boltzmann relaxation parameter is set tor = 1, and thus 
the time step scales like At oc Aa; 2 (diffusive scaling). 
The simulation box is a cube with H 3 fluid lattice nodes. 
The capsule with initial radius r is placed at the center of 
the box. We choose the x-axis as the velocity direction 
and the j/-axis as the vorticity direction. The velocity 
gradient is along the z-axis. The bounding plates are 
located at z = ±H/2. 

In Sec. IIV A[ we first determine the effects of the rel- 
ative simulation box size H/r, the Reynolds number Re, 
the reduced shear rate G, and the LBM relaxation pa- 
rameter r on the deformation state of the capsule. Based 
on those first results, the influence of the membrane tes- 
sellation is analyzed in Sec. IIVB1 In Sec. IIV CI we test 
the influence of the interpolation stencils on the numer- 
ical results and specify the convergence behavior of the 
IBLBFEM scheme. 



A. General simulation parameters 

We note that, in the limit Re ^ 1, H/r ^> 1, and 
G -C 1, all curves D(nt)/G, with k = j/G, collapse on 
the same master curve. This can be used to study the 
impact of finite size effects (H/r ^> 1), inertia (Re 1), 
nonlinear contributions (G 5^ 1) from the constitutive 
membrane model, and the LBM relaxation parameter t 
by inspecting the deviations from the master curve. For 
all simulations in this section, the 4-point stencil ^4, Eq. 
(fl"8|l. and an icosahedron-based mesh with Nf — 1280 
faces and r = 5 have been used. 

Similar to Sui ct al. (l2j . we have first determined 
the minimum system size in terms of H/r to safely ne- 
glect self-interaction of the capsule or interactions with 
the walls. The simulation parameters are Re = 0.02, 
G = 0.01, and r = 1. We have examined the cases 
H/r = 6,8,10, and 12. The resulting time evolution 
of the Taylor parameter is shown in Fi 
to the same conclusion as Sui et al 
of H/r = 10 is sufficient for modeling unbounded simple 



7(a) 



We come 
that a box size 



shear flow. The difference between the plateau values of 
the Taylor parameter for H/r = 10 and 12 is less than 
0.5%. For this reason, the system size is taken to be 
H/r = 10 in all following simulations. However, it is 
obvious that the deformation parameters are too large. 
The expected value, D/G = 6.25, is shown in Fig. [71 We 
will see in this section that this is not related to effects 
caused by inertia, nonlinear membrane response, or the 



relaxation parameter. This observation will be discussed 
in more detail in the following sections, and the explana- 
tion will be given in Sec. IIV CI 

In the next step, we have investigated the validity of 
the Stokes flow assumption. A series of simulations with 
H/r = 10, G = 0.01, and r = 1 has been carried out. The 
tested Reynolds numbers arc Re = 0.01, 0.02, 0.05, and 

0. 1. Sui et al. [lj| have observed that inertia effects are 
small up to a Reynolds number of about 0.025. We come 
to a similar conclusion. In Fig. |7(b)[ the deformation 
parameter is shown for different Re. Differences between 
Re = 0.01 and 0.02 can only be noticed slightly in the 
transient evolution. Since we are interested in the steady- 
state behavior, the Reynolds number will be set to 0.02 
in all subsequent simulations. We have also studied the 
time evolution if the equilibrium in Eq. @ is linearized, 

1. e., substituted by 



/,- q = Wip(l + 3c 4 • u) 



(31) 



One can show that this leads to the removal of the ad- 
vection term u ■ Vu in the Navier-Stokes equations. Still, 
the partial time derivative du/dt is present. Within the 
kinematic framework of LBM, the time derivative cannot 
be removed, and exact Stokes flow cannot be simulated. 
We have observed that the time evolution of the defor- 
mation parameter does not change when the advection 
term is removed with respect to the cases where the full 
equilibrium has been considered (data not shown). We 
thus assume that the Reynolds number effects visible in 
Fig. |7(b)| are due to the partial time derivative and not 
the advective term. Still, neglecting the second-order 
term in the equilibrium, the computing time for LBM 
could be decreased significantly (« 25%). This observa- 
tion can be useful for high performance simulations at 
small Reynolds numbers. However, since we have tested 
the first-order equilibrium afterwards, for the remaining 
simulations the second-order equilibrium has been em- 
ployed. 

Additionally, we have tested up to which value of the 
reduced shear rate G the linearity assumption, Eq. (f26|) . 
is valid. In Fig. 7(c) the time evolution of D is shown. 
The simulation parameters are Re = 0.02, H/r = 10, 
and t = 1. The reduced shear rates arc G = 0.005, 
0.01, 0.02, and 0.04. We find that G = 0.01 is sufficient 
to ensure the validity of the linear approximation. For 
G > 0.02, the deviations become significant, and second- 
order correction terms should be included. In all the 
following simulation, we have chosen G = 0.01. Note 
that the plateau values of D/G between G = 0.005 and 
0.01 differ by less than 2%. 

It is known that the BGK LBM relaxation parameter 
r plays a critical role in the correct setup of the simula- 
tions. Both the accuracy of the bulk LBM and the no- 
slip bounce-back boundary conditions depend on its value 
13 13 HI! • Recently, Le and Zhang [43] reported that 
combined IBM-LBM simulations are strongly affected by 
the magnitude of r, especially if r > 1. We address this 
issue by comparing simulations with different values of 
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Figure 7. Influence of finite size, inertia, nonlinear membrane response, and LBM relaxation parameter on the deformation 
parameter, which is shown as a function of the reduced time Kt, where k = j/G. |(a)| finite size: Re = 0.02, G = 0.01, r = 1, 
and H/r = 6, 8, 10, 12; [(b)] inertia: H/r = 10, G = 0.01, r = 1, and Re = 0.01, 0.02, 0.05, 0.1; \(cj\ nonlinear membrane response: 
Re = 0.02, H/r = 10, r = 1, and G = 0.005,0.01,0.02,0.04; [(d)] relaxation parameter: Re = 0.02, H/r = 10, G = 0.01, 
and t = 0.6,0.8, 1.0, 1.5,2.0. In all simulations, Nf = 1280, r = 5, and <f>4. have been used. The influence of the finite box is 
negligible for H/r > 10, and inertia effects are unimportant for Re < 0.02. The inset in subfigure | (c) | shows the final value of 
D/G versus the reduced shear rate G. The dashed line is the analytic result for small deformations. For G > 0.02, nonlinear 
effects become obvious, but G = 0.01 is a good approximation for a small deformation. The relaxation parameter r should 
not be much larger than unity. It can be seen from the inset of subfigure (d) that the final value of D/G strongly depends 
on r, if r > 1. The theoretical value, D/G = 6.25, is shown in all subfigures. One notices that the deformation parameters 
consequently are too large. 



the relaxation parameter for Re = 0.02, H/r = 10, and 
G = 0.01. The results are shown in Fig. |7(d)| It can 
clearly be seen that the results are relatively indepen- 
dent of t for r < 1. If r becomes significantly larger than 
unity, the solutions start to diverge. This effect may be 
related to the single relaxation time (BGK) LBM which 
has been employed. Having the efficiency of the simula- 
tions in mind, it is desired to keep the simulation time 
as short as possible. This can be achieved by increasing 
the time step At. From 

r - 1/2 Ax 2 , . 

v = " T- (32) 

3 At ' 

we can see that increasing r also increases At, if the 



kinematic viscosity v and the spatial resolution Ax are 
kept fixed. This reduces the number of necessary time 
steps in the simulations. However, if r becomes too large, 
drastic numerical artifacts appear, and the simulations 
become unreliable. It seems to be a compromise to choose 
r 1. 

Summing up, we come to the first conclusion that a 
cubic box of length H/r = 10, a Reynolds number of 
Rc = 0.02, a reduced shear rate of G = 0.01, and a LBM 
relaxation parameter of r = 1 are excellent approxima- 
tions to the unbound Couette flow at vanishing Reynolds 
number in the linear elastic limit. However, the deforma- 
tion parameters are larger than expected. Theory pre- 
dicts D/G = 25/4 = 6.25 for the Skalak membrane with 
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k s = k a , but typical values are D/G ~ 7.0 in the pre- 
sented simulations. The expectation value is also shown 
in Fig. [7] There are two possible reasons for this dis- 
crepancy: first, the membrane mesh may be too coarse. 
Second, the interpolation and spreading, required for the 
coupling between Lagrangian and Eulcrian meshes, may 
introduce numerical artifacts which effectively lead to a 
softer or, equivalcntly, larger capsule. We stress that the 
IBM generates an artificial length scale Li related to the 
width of the interpolation. If this length scale is not small 
with respect to the length Lm of significant changes in 
fluid velocity and membrane tensions, a detrimental ef- 
fect of the interpolation and spreading is expected. In 
Sec. IIV B\ we first test the effect of the mesh discretiza- 
tion. Afterwards, in Sees. IIV C\ we turn to the influence 
of the IBM interpolation stencil on the simulations. 



Mesh discretization 



In Sec. IIV Al we have stated that the deformation pa- 
rameters are larger than expected from linear theory. 
From the discussions there, finite size effects of the sim- 
ulation box, nonlinear membrane responses, and inertia 
effects could be excluded. Also the choice of the LBM 
relaxation parameter cannot be responsible for the devi- 
ations. The most probable explanations are spatial dis- 
cretizations due to the membrane tessellation (cf. Sec. 
Ill D|) or the discrete delta functions for interpolation and 
spreading, Eqs. p6|) -(fT8 |) . Since length and time scales in 
the LBM are strongly coupled 0, Eg], spatial and tempo- 
ral discretization errors cannot be studied independently. 

In order to test the influence of the details of the mesh 
resolution, we have performed three simulations with 
identical parameters (Re = 0.02, H/r = 10, G = 0.01, 
r = 1, r = 5, and 04 ), but different meshes (icosahedron- 
based, N f = 320, 1280 and 5120 faces). As a conse- 
quence, the average distances between neighboring mesh 
nodes (i.e., the average edge length of the face elements) 
differ: 1/ Ax = 1.50 for the coarsest, 0.75 for the interme- 
diate, and 0.38 for the finest mesh. The results for the 
deformation parameter, the inclination angle, and the an- 
gular velocity are shown in Fig. [8j All curves for D nearly 
collapse. Thus, the mesh discretization can be dropped 
as explanation for the numerical softening of the capsules 
since the difference between the meshes with Nf = 320 
and 5120 faces is small. The results for the inclination 
angles are very similar, but all of them show a deviation 
from the expected value as well. We will come back to 
this point in Sec. IIV CI A good approximation of the an- 
gular velocity uj can only be provided by a large number 
of face elements. 

At this point, we do not see any reason why the av- 
erage distance between nodes should be less than Ax/2, 
which is sometimes claimed in the literature (26|. This 
result is very important from an efficiency point of view. 
In a simulation of a dense suspension of dcformable parti- 
cles, most of the computing time is required for the IBM 



interpolation and spreading. Since the number of IBM 
calculations is proportional to the number of mesh nodes, 
it is worth examining to which extent the membrane res- 
olution can be reduced at fixed fluid lattice resolution 
without significantly decreasing the accuracy of the sim- 
ulations. As a consequence, the computing time could be 
decreased. It has to be noted that, independent of this 
result, a finer mesh has to be used when the radius is sig- 
nificantly increased since the average distance between 
neighboring nodes would become too large eventually. 

Additional tests have been performed to analyze the 
effect of different mesh tessellations (cf. Sec. IIID|) . The 
simulation parameters are as above, but Nf = 1280 
(icosahedron-based), 1284 (Gmsh), and 1278 (CGAL). 
The results are collected in Fig. [9l There is no difference 
in the evolution of the deformation parameter, but the 
inclination angle is not correctly captured at small de- 
formations when the Gmsh- and CGAL meshes arc used. 
Strong deviations can be seen in the time evolution of 
the angular velocity for the mesh created by Gmsh. The 
reason is that the computation of uj is most susceptible 
to numerical artifacts caused by reduced homogeneity of 
the mesh. However, the general behavior of the capsules 
is similar, and the deformation state is not strongly in- 
fluenced by the details of the mesh tessellation. 

Concluding this section, it is found that details of the 
mesh (tessellation method and resolution) do not sig- 
nificantly change the deformation state of the capsule. 
Thus, the deviation of the deformation parameter from 
the expected analytic value cannot be caused by the dis- 
cretization of the capsule membrane. Even a small mesh 
resolution is sufficient to accurately describe its deforma- 
tion behavior, at least at small values of G, cf. Fig. 8(a)| 
Some details of the capsule are not correctly captured 
by the Gmsh and CGAL meshes. We will employ the 
icosahedron-based mesh in all remaining simulations. 



C. Interpolation and spreading 

We have examined the influence of the discrete mem- 
brane mesh on the simulations in Sec. IIV Bl where only the 
4-point interpolation stencil 4>i has been employed. The 
conclusion is that the discrepancy between the observed 
and predicted steady-state values of the deformation pa- 
rameter D and inclination angle should be related to 
the interpolation stencil of the IBM. This open point will 
now be discussed in more detail. For convenience and 
clarity, we define a numerical capsule radius N r = r/ Ax, 
indicating the number of fluid lattice nodes covered by 
the actual radius r. 

A consequence of the presence of the interpolation 
stencils, Eqs. (fl6)) (fT8|) . is the finite numerical width 
of the membrane. Consequently, the computed solu- 
tions should converge to the analytic predictions for 
Lj/Lm 0. Here, Lj is the length scale associated 
with the numerical membrane thickness, and Lm can be 
regarded as the radius of the membrane. The values of Lj 
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Figure 8. Effect of the icosahedron-based mesh resolution on (a) the deformation parameter, (b) the inclination angle, and 
|(c) | the angular velocity. The simulation parameters are Re = 0.02, H/r = 10, G = 0.01, r = 1, and r = 5. The 4-point 
interpolation stencil has been used. There is virtually no change in the deformation parameter when the 5120-element mesh is 
replaced by the 1280- or the 320-element mesh. The same is valid for the inclination angle. The angular velocity is extremely 
sensitive to the number of face elements, but convergent behavior is evident. The expected analytic values for the deformation 
parameter and the inclination angle are also shown in the subfigures. 



for the available interpolation stencils may be identified 
with Ace, 3Acc/2, and 2 Ax for 02, (f>3, and 04, respec- 
tively. Since numerical width effects can also be seen in 
multiphase simulations, we refer to a review about diffuse 
interface methods by Anderson et al. [fij) • 

As discussed before, one can find comments on the ra- 
tio of mesh and lattice resolution 1/ Ax in the literature. 
Peskin (2(| suggests that two adjacent nodes on the mem- 
brane mesh should have a mean distance / < Ax/ 2. This 
way, it is argued, no fluid could leak between the 'holes' in 
the mesh. However, a more stringent motivation is miss- 
ing. In Sec. lIVBl we have not found evidence supporting 
this claim. The results indicate that also an average node 
distance of I /Ax = 1.5 may lead to a reasonable accu- 
racy. We stress that a necessity for smaller distances I 
may arise in the limit of stiff particles or strongly de- 
formed membranes. 

Beside the choice of the interpolation stencil 0, the ra- 
tio 1/ Ax of the mesh and lattice resolutions is the only 



freedom left in the IBM. If 1/ Ax becomes too large, fluid 
will eventually leak through the membrane. On the other 
hand, if / is chosen to be too small, the spatial hydrody- 
namic resolution becomes worse since less fluid lattice 
nodes cover the capsule, and the ratio Li/Lm increases. 
Within the framework of IBM, there is no obvious, sim- 
ple way of estimating the optimum value of 1/ Ax. In- 
tuitively, both limiting cases restrict the accuracy of the 
simulations. Therefore, it is important to understand in 
which range it is safe to operate. 

In order to quantify the dependence of the deforma- 
tion parameter D and the inclination angle 9 on Lj/ Lm , 
we have performed two different kinds of studies, each 
employing 2 , 03, and 04 . The simulation parameters 
always are Re = 0.02, G = 0.01, H/r = 10, and r = 1. 

a. Convergence for fixed mesh resolution Nf In this 
simulation series, we study the influence of the hydrody- 
namic resolution N r alone, i.e., we keep the mesh reso- 
lution Nf constant, and 1/ Ax changes. This way, it is 
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Figure 9. Effect of the tessellation method on (a) the deformation parameter, (b) the inclination angle, and (c) the angular 
velocity. The simulation parameters are the same as those in Fig. [8] but Nf = 1280 (icosahedron-based), 1284 (Gmsh), and 
1278 (CGAL). The deformation parameter is not influenced by the type of the mesh, but the inclination angle and the angular 
velocity are detrimentally affected. The inclination angle is shown at early times only to reveal the deviations. Due to the 
reduced homogeneity of the meshes, the correct inclination angle is not captured at small deformations where it should be 
9/tt « 0.25. Also the angular velocity shows unphysical behavior, especially for the mesh created with Gmsh. 



possible to study the effect of a varying ratio I /Ax. The 
employed mesh resolutions are Nf = 320 and 1280. 

For the mesh with 320 faces, we have tested N r = 3, 
4, 5, and 6, corresponding to l/Ax — 0.90, 1.20, 1.50, 
and 1.80, respectively. The results are shown in Fig. [TO] 
Although the mesh resolution Nf is small and the mesh 
ratio 1/ Ax significantly larger than 0.5, it can be seen 
that the physics of the system is roughly captured. The 
accuracy of the solutions increases with a larger magni- 
tude of N r . For <p2, fluctuations are visible. The smallest 
radius, N r = 3, yields inaccurate results. In this case, 
Li and Lm are comparable. It is interesting to note that 
even for 1/ Ax = 1.80 no detrimental effects appear. The 
fluid still does not seem to penetrate the membrane. 

Additionally, we have tested the mesh with 1280 faces 
and N r = 3, 5, 7, and 9, corresponding to 1/ Ax = 0.45, 
0.75, 1.06, and 1.36, respectively. The results arc shown 
in Fig. 111! The overall behavior of the numerical results is 
similar to those of the series with Nf — 320. The smallest 



radius, N r = 3, is less accurate, and no penetration of the 
fluid is visible in the presented parameter range. 

We have observed that the 2-point interpolation sten- 
cil 4>2 fails when I /Ax > 2 (data not shown). At this 
point, the spacing between neighboring mesh nodes is 
so large that fluid can penetrate the capsule membrane. 
For </>3 and (j>± , we have not observed a similar behavior at 
1/ Ax = 2. The probable explanation is the larger range 
of the interpolations, still keeping the fluid from passing 
through the membrane. 

The above studies strongly suggest that the mesh ra- 
tio can be safely chosen somewhere between 0.5 and 1.5 
without compromising the impermeability of the capsule. 
This is an important result since it allows us to reduce the 
computational requirements by a proper choice of 1/ Ax 
without loss of accuracy. 

b. Convergence for fixed mesh ratio I /Ax In this 
second series, we investigate the coupled convergence 
when both the mesh and the hydrodynamic resolutions 
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Figure 10. Behavior for varying mesh ratios I /Ax for a fixed mesh with Nf = 320 faces with the interpolation stencils (j>2 
in subfigures (a) and (b) (f>3 m subfigures (c) and (d) and (f>4 m subfigures (e) and (f) N r = 3, 4, 5, and 6 correspond to 



l/Ax = 0.90, 1.20, 1.50, and 1.80, respectively. 



are increased by the same rate, i.e., the mesh ratio I /Ax 
is fixed. The mesh and hydrodynamic resolutions are 
Nf = 1280 and H = 35, N f = 5120 and H = 70, and 
Nf = 20480 and H = 140, respectively. Note that for 



Nf = 320, a mesh ratio of 1/ Ax = 0.53 leads to quite in- 
acceptable results. This is closely related to the fact that 
the hydrodynamic radius becomes too small compared 
to the numerical width of the membrane. Here, we see 
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Figure 11. Behavior for varying mesh ratios I /Ax for a fixed mesh with Nf = 1280 

and 



in subfigures (a) and (b) 
l/Ax 



in subfigures (c) and (d) 



0.45, 0.75, 1.06, and 1.36, respectively. 



hi in subfigures (e) and (f) 



faces with the interpolation stencils (j>2 
N r ~ 3, 5, 7, and 9 correspond to 



again that using the freedom in the choice of I /Ax allows 
us to use mesh resolutions which would not be available 
otherwise. The mesh ratio is I /Ax = 0.53 in all cases. 
The results are shown in Fig. [T2J It is obvious that the 



steady-state magnitudes of D and 9 converge to their an- 
alytic values (D a /G = 6.25 and a /ir ps 0.231). In order 
to quantify the results, the errors at Kb = 120 are listed 
in Tab. HI1 and graphically shown in Fig. [T3] The conver- 
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Convergence of the IBLBFEM for the interpolation stencils (f>2 in subfigures | (a) | and | (b) | 4>s in subfigures | (c) | and | (d) | 
The mesh ratio 1/ Ax is 0.53 in all simulations. The mesh and fluid resolutions are Nf = 1280 
A?7= 5120 and H = 70 (dashed lines), and Nf = 20480 and H = 140 (solid lines), respectively. 
Convergence to the analytic predictions is observed in all cases. The numerical errors at nt = 120 (ft = j/G) are also shown in 
Tab. HI] and Fig. [H 



Figure 12 
and 4>4 in subfigures (e) and (f) 
and H = 35 (dotted lines) 
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Table II. Convergence of the IBLBFEM for the interpolation stencils 02, 03, and 04. The deviations of the deformation 
parameter D and the inclination angle 9 at nt — yt/G = 120 are shown. The relative deviations are defined as SD = 
(D a — Da) /D a and 89° = (9° — 0°)/0 a . The subscripts s and a denote 'simulation' and 'analytic'. The angle 8° oc G is the 
opposite angle to 9, defined in Eq. (|28[) . The convergence order a is taken from a fit to the function 8D, 59° oc N~ a . For 88° 
and 02, a meaningful convergence order could not be obtained due to early mesh degradation. A graphic representation of this 
table is shown in Fig. 1131 
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Figure 13. Convergence of the IBLBFEM for the interpolation stencils 02, 03, and 04. In subfigure (a) the error of the 
deformation parameter, 8D, is shown for increasing mesh and fluid resolutions with fixed 1/ Ax = 0.53. The analog results for 
the error of the inclination angle, 89°, are shown in subfigure (b) The data is taken from Tab. ILT1 



gence order is clearly better than 1, cf. Tab. HT1 The only 
exception is the inclination angle with 02 ■ This is caused 
by mesh degradation, cf. Sec. IIVDI Since the LBM is 
second-order accurate and the IBM for sharp interfaces 
formally first order, the convergence order of the coupled 
system should be between 1 and 2 as observed here. Un- 
fortunately, we are not aware of any theory which could 
predict the convergence behavior of the coupled system. 

c. Effective deformability and rescaling From the re- 
sults in this section, we draw two main observations. 

First, the numerical magnitudes of the deformation pa- 
rameter D and the inclination angle 9 approach the ex- 
pected values when the hydrodynamic resolution N r is 
increased. Within the valid region of I /Ax, this state- 
ment also holds if only the hydrodynamic resolution is 
increased and the mesh resolution is kept constant. How- 
ever, the formally correct approach is to gradually refine 
both meshes simultaneously. The IBLBFEM accurately 
captures the physics of dcformablc capsules in an ambient 
fluid within the chosen parameter ranges and in the limit 
of infinite resolution. The convergence order is between 



1.5 and 2, cf. Tab. [0] and Fig. ED 

Second, we conclude that the deviations between ex- 
pected and analytical values are an IBM artifact. The 
average deviations are usually smaller for a narrower in- 
terpolation stencil, cf. Tab. [TT| (except 02 at large resolu- 
tion). This supports the idea that the numerical width 
of the interpolation stencil affects the deformation be- 
havior of the capsule. The reason for the effect of 02 at 
large resolutions is an accelerated mesh degradation. We 
will come back to this point in Sec. IIVDI On the first 
glance, a narrower interpolation stencil does a better job 
in capturing the physics of the problem. However, one 
finds that fluctuations are more pronounced when 03 and 
especially 02 are employed since those interpolations are 
not as smooth as 04. This can be seen in the plots of 9 in 
Figs. [TU1 HU and 1121 Making a good choice for the inter- 
polation stencil means balancing the strengths and weak- 
nesses of those stencils. We have also seen that, even for 
coarser resolutions, the numerical results are sound and 
reliable, as long as the presence of the finite membrane 
thickness is properly taken into account. 
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Table III. Effective reduced shear rates G/G and effective errors 56*° = (8° — 6°(G))/9°(G) of the opposite angle. G is defined 
in such a way that the deformation parameter has no error, i.e., D a = D a (G). The effective deviations 56° are clearly reduced 
compared to the raw data in Tab. [TTJ 
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It is worthwhile to discuss the deviations of the defor- 
mation parameter D and the opposite inclination angle 
9°, Eq. [251 in more details. While the deviations SD and 
89 become smaller when the hydrodynamic resolution N r 
is increased, SD and 59° are always of the same order. 
Since D oc G and 9° oc G, this enables us to introduce an 
effective reduced shear rate G to partially counteract the 
effect of the finite membrane thickness due to the pres- 
ence of the interpolation stencil. Starting from Eq. (|23)l . 
the effective reduced shear rate can be written as 

G=?f, (33) 

i.e., the transition G — > G can be due to a transition 
r — > r or k s — > k s . We claim that the fluid properties 
7 and rj are well defined, thus, we do not allow their ef- 
fective modification. In the present simulations, we have 
the choice to define an effective radius f or an effective 
stiffness k s in such a way that the simulation data are 
more accurate. It is left for future research whether a 
redefinition of r or k s is more useful. In the present 
simulations, both approaches are equivalent, but due to 
particle-particle interactions in simulations with multi- 
ple capsules, the effect of the finite membrane width is 
more complex and a simple redefinition of r or k s may be 
not straightforward. In general, the rescaling factor G/G 
is a function of the interpolation stencil and the hy- 
drodynamic and the mesh resolutions, N r and Nf. The 
factor G/G can be obtained from Tab. [TTJ In order to 
investigate this idea further, we have defined G for each 
simulation in Tab. [Tj] in such a way that the deviation 
SD := (D s - D a (G))/D a (G) is identically zero. The cor- 
responding error 89° := (9° - 9°(G))/9°(G) is shown in 
Tab. IIIII All effective errors S9° are considerably smaller 
than the original errors given in Tab. [TTJ stating that the 
redefinition of G can be used to compensate — at least 
partially — the numerical width effect due to the interpo- 
lation stencil. It is expected that for interacting capsules 
with more complex shape, e.g., red blood cells, a sim- 
ple rescaling is not as straightforward as for an isolated 
spherical capsule. 

To sum up the above findings, we note that the domi- 
nant source of numerical deviation from the analytic so- 
lution in the combined IBLBFEM is the presence of the 
interpolation stencils. They introduce a numerical width 
of the membrane, leading to slightly modified physical 



behavior. Those deviations could in principle be counter- 
acted by introducing an effective radius or effective stiff- 
ness of the capsule. In the limit of infinite resolution, the 
IBLFFEM accurately captures the physics of the coupled 
system of capsule and fluid. The influence of the choice 
of l/Ax is small over a wide range, 0.5 < l/Ax < 1.5. 
This introduces a freedom to the numerical implemen- 
tation which can be used to decrease the computational 
cost of the simulations. 



D. Mesh degradation and volume drift 

In this section, we will address some issues which have 
not been covered in the discussion about the mesh influ- 
ence and the convergence studies in Sees. IIV AIITVCl 

Although we have claimed that the mesh resolution 
does not have a significant influence on the results (cf. 
Sec. IIV Bp . we have reported different deviations of the 
deformation parameter D and the inclination angle 9 
when 4>2 and varying mesh resolutions are used (cf. Tab. 
ITU) . The reason is that the values in Tab. [TT] are taken 
at time nt = 120 where the mesh has already started to 
degrade when 02 is employed. This effect seems to be 
most severe when the hydrodynamic resolution is small- 
est (i.e., small radius N r and large number of faces Nf). 
A similar behavior is only weakly noticeable for 03 and 
04, indicating that those interpolation stencils do a bet- 
ter job in preserving the mesh. This observation is an 
indication that the average node distance l/Ax should 
not be too small. The mesh degradation can be captured 
by computing the average face element area and its stan- 
dard deviation. For a degrading mesh, one expects those 
values to leave the steady state gradually. Our findings 
are illustrated in Fig. [Ml Only a few curves are shown. 
The major observation is that increasing the mesh resolu- 
tion but keeping the hydrodynamic resolution unchanged 
leads to an accelerated and undesired mesh degradation. 
Although the physical behavior of the capsule is indepen- 
dent of the mesh resolution at small times (cf. Sec. IIV Bl) . 
numerical artifacts become progressively more important 
at later times. For long-time simulations, the shear en- 
ergy W may be not sufficient to control the mesh. 

Taking the results of the simulations presented in Sec. 
IIV CI we have observed typical volume deviations SV/Vq 
between 2 • 10 -4 and 8 • 10 -4 for 02, between 5 • 10 -5 and 
2 ■ 10~ 4 for 3 , and between 9 ■ 10~ 6 and 3 • 10" 5 for 4 



20 




Figure 14. Time evolution of average and standard deviation of the relative face area deviation 5A/A. At t = 0, all faces 
are undeformed and hence 8A/A = and ctsa/a = 0. Due to the deformation of the capsule, (a) the average area deviation 
SA/A increases until it reaches a steady state. The steady state eventually collapses when the mesh degrades. This is most 
significant for </>2 at the highest resolution of the membrane mesh (Nf = 20480). The degradation is much less significant 
for cj>2 and Nf = 5120, indicating that too fine a mesh (i.e., too coarse a hydrodynamic resolution) is detrimental. |(b)| The 
standard deviation asA/A of the relative area deviation increases with time, especially for a mesh with high resolution but poor 
hydrodynamic resolution. 



at nt = 120. As expected, the largest deviations corre- 
spond to the smallest radius N r and vice versa. Clearly, 
the 4-point interpolation function is much more appro- 
priate to control the capsule volume without taking addi- 
tional measures. A higher hydrodynamic resolution also 
reduces the volume drift . 

In the present simulations, the volume drift is negli- 
gible. However, in long-time simulations with smaller 
resolution, the volume drift could become a significant 
problem. Taking into account restoring forces originat- 
ing from volume, surface, and bending energies may help 
to avoid this deficiency. 



V. CONCLUSIONS 

We found that the choice of the LBM relaxation pa- 
rameter r can have a detrimental effect on the combined 
IBM-LBM simulations. As long as r < 1, the simula- 
tion results barely depend on the actual magnitude of r. 
For larger relaxation parameters, strong deviations can 
be observed, and the simulation results obtained by dif- 
ferent r are no longer comparable. On the other hand, 
we have seen that the choice of r strongly affects the 
simulation time since r can be used to change the LBM 
time step. We have used r = 1 in the simulations as a 
compromise between accuracy and efficiency. 

The common assumption that the average node dis- 
tance I should be smaller than half a fluid lattice con- 
stant Ax could not be supported in our investigations. 
This leads to the important consequence that, at least 
within a certain range, the hydrodynamic resolution and 
the mesh resolution can be changed independently with- 
out significantly changing the physics in the simulations. 



Due to its locality and purely algebraic structure, the 
LBM algorithm is fast and efficient. It is well known that 
the computing speed of pure LBM simulations is usually 
restricted by memory access, but not by CPU power. 
Considering the capsule immersed in the fluid, the com- 
putation of the membrane forces, the velocity interpola- 
tion, and the force spreading is computationally more de- 
manding and basically limited by the CPU power. In the 
present simulations, the fluid volume is large compared to 
the membrane area, and nearly all the computing time 
is consumed by the LBM component. For simulations 
with moderate or high particle volume fractions, how- 
ever, the IBM consumes most of the computer resources. 
We have observed that the computing time for the ve- 
locity interpolation and spreading of the forces is much 
larger than that for the computation of the membrane 
forces itself. The number of IBM interpolations at each 
time step is 2m d N n in d dimensions, using a stencil <p m 
with a support of m lattice nodes along each direction 
and a total of N n membrane nodes. One way to save 
computing time is to choose an interpolation stencil with 
smaller support. Although the relevant physics seems 
to be captured with the 2-point stencil as well, stronger 
fluctuations are introduced, and the mesh undergoes an 
accelerated degradation (cf. Sees. HVCl and UVDp . How- 
ever, a smaller support is also equivalent to a smaller 
numerical thickness of the membranes, which could be of 
advantage in a dense suspension of dcformable particles 
where the distances between the membranes are small. 

The major discretization error is introduced by the 
IBM interpolation stencil <p, but not by the discretization 
of the membrane itself. The hydrodynamic resolution is 
more important than the mesh resolution. This observa- 
tion leads to the following conclusion: if the accuracy of 
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the simulation shall be increased, it is the optimum ap- 
proach to first increase the number of fluid lattice nodes 
for fixed mesh resolution. If on the other hand the com- 
puting time of the simulation should be decreased, the 
hydrodynamic resolution should be kept and the mesh 
resolution be reduced. 

We observed that the presence of the interpolation 
stencils <f> effectively changes the membrane properties of 
the capsules due to the numerical thickness of the mem- 
brane. This can be captured by defining an effective re- 
duced shear rate G by considering an effective capsule 
radius f or an effective stiffness k s . The numerical mem- 
brane width decreases when the hydrodynamic resolution 
is refined, i.e., when the number N r of fluid lattice nodes 
covered by the capsule radius is increased. The numer- 
ical results converge to the analytically expected values 
for large N r . Consequently, making use of the knowledge 
of an effective radius r or stiffness k s in principle allows 
for the reduction of the hydrodynamic resolution with- 
out a significant loss of accuracy. This approach may be 
used to massively save computing time in simulations of 
multiple deformable particles immersed in a fluid. 

We stress that, in this paper, only small deformations 
of a single capsule have been considered. The mesh res- 
olution plays only a minor role in this case. However, 



the number of mesh points may be of higher importance 
in simulations with large membrane deformations since 
the faces remain always flat. For that reason, the analy- 
sis of the influence of the mesh resolution should also be 
performed for large deformations. The effect of particle- 
particle interactions in dense suspensions on an effective 
radius or stiffness caused by the interpolation stencils is 
also not obvious at this point. Those investigations are 
left for future research. 

The present paper contains new contributions: regard- 
ing the effect of the interpolation stencils on the capsules' 
deformation behavior, the significance of the hydrody- 
namic resolution compared to the mesh resolution, and 
how those insights may be used to boost the efficiency 
of the related simulations. These results are hoped to be 
useful for the simulation of multiple deformable objects 
immersed in a fluid. 
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